#include "ParallelJacobisMethod.h"

//#define PJM_DEBUG_FCK 
//#define PJM_DEBUG_JM 
//#define PJM_DEBUG_INPUT 
//#define PJM_DEBUG_OUTPUT 

int findMyChunk(
               struct extendedMatrix *M, /* IN */
               struct unknownsStrip *X,  /* OUT */
               int rank,                 /* INT */
               int nProcess              /* INT */
               ){

    unsigned int D, L, U;

    /*
    Exemple: 10 lines 4 process;
    Process0 is the master, so whe skip it

    Process1 takes lines 0, 1, 2
    Process2 takes lines 3, 4, 5 
    Process3 takes lines 6, 7, 8, 9
    
    D = 10/3 = 3.33
    
    L - is multiple of D
    U - Next rank minus 1

    if the process is the last

    L - is multiple of D
    U - last Line
    */
    
    D = (M->J_ORDER)/(nProcess-1); 
    L = D*(rank-1);

    #ifdef PJM_DEBUG_FCK 
        printf("RANK: %d NPROCESS %d \n", (rank-1), nProcess);
    #endif

    if(rank!=(nProcess-1)) U = (D*(rank))-1;
    else U = (M->J_ORDER-1);
     
    X->lowerChunk = L;
    X->upperChunk = U;

    #ifdef PJM_DEBUG_FCK
        printf("CHUNK: L=%u U=%u  U-L=%u D=%u\n", L, U, U-L, D);
    #endif

    return (U-L);
}

struct unknownsStrip *ParallelJacobisMethod(
                    struct extendedMatrix *M,
                    struct unknownsStrip *X
                    ){

    unsigned long i;
    unsigned long L, U, N; 
    struct unknownsStrip *XN;    
    
    XN = initStrip();

    XN->lowerChunk = X->lowerChunk;
    XN->upperChunk = X->upperChunk;

    L = X->lowerChunk;
    U = X->upperChunk;
    N = M->J_ORDER;

    XN->length = (U-L)+1;
    XN->vector = (double*) calloc(XN->length, sizeof(double));
    XN->isItok = 1;

    #ifdef PJM_DEBUG_INPUT
        printf("\n\nX = %lu %lu  XN->length= %u\n", L, U, XN->length);
        for(i=0; i<X->length; i++)
            printf("%.3lf ", X->vector[i]);
        printf("\n\nB = %lu %lu\n", L, U);
        for(i=0; i<X->length; i++)
            printf("%02.3lf ", M->vector[i]);
        printf("\n");
    #endif

    for(i=0; i<XN->length; i++){
        XN->vector[i] = jmFindXi(M, i+L, X->vector);
    }

    #ifdef PJM_DEBUG_OUTPUT
        printf("\n\nXN = %lu %lu\n", L, U);
        for(i=0; i<X->length; i++)
            printf("%.3lf ", XN->vector[i]);
    #endif

    return XN;
}

double pjmFindXi(struct extendedMatrix *M,
                  unsigned long line,
                  double *X){
    double Xi, SUM=0, Ad, Aij;
    unsigned long j;

    // Xi = b
    Xi = M->vector[line];    

// TODO
    #pragma omp parallel for reduction (+:SUM)
    for(j=0; j<M->J_ORDER; j++){
        Aij = M->matrix[line][j];
        if(j==line){
            // Element from the diagonal
            Ad = Aij;
            continue;
        }else{
            // Not from the diagonal
            // Xi = b-SUM(Aij*Xk)
            SUM += Aij*X[j];
        }
    }
    
    Xi -= SUM;
    Xi/=Ad;
            /*  Xk+1 = b-SUM(Aij*Xk)
                       -------------
                           Aii        */
    return Xi;
}

struct unknownsStrip *initStrip(){

    struct unknownsStrip *X;
    X=(struct unknownsStrip*) malloc(sizeof(struct  unknownsStrip));

    X->length=0;
    X->lowerChunk=0;
    X->upperChunk=0;
    X->isItok=0;
    X->vector=NULL;

    return X;
}

void freeStrip(struct unknownsStrip *X){

    if(X==NULL) return;
    if(!X->isItok){
         free(X);
         return;
    }
    free(X->vector);
    free(X);
}
